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ABSTRACT 
The  partial  differential  equation  for  heat  diffusion  is  numerically 
integrated  by  the  Runge-Kutta-Cill  method.  A  solution  is  obtained  for 
the  diurnal  temperature  variation  with  a  hounded  coefficient  of  eddy 
diffusivity  vhich  varies  periodically  with  time  and  nonlinearly 
with  height.  The  surface  wave  is  represented  by  the  sum  of  a 
diurnal  and  a  semidiurnal  harmonic  wave.  The  results  may  be  interpreted 
to  apply  over  a  fairly  broad  range  of  diffusivity  values  and  height. 
With  appropriate  choices  of  the  various  parameters,  reasonably  good 
agreement  is  obtained  between  theoretical  and  observational  values  of 
amplitude  reduction  and  phase  lag. 

The  author  wishes  to  express  his  appreciation  for  the  assistance 
and  encouragement  of  Professors  G.  J.  Haltiner  and  E.  J.  Stewart  of 
the  U.  S.  Naval  Postgraduate  School  in  this  investigation. 
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1.   Introduct Ion. 

Of     --y,   problems  facj      ^orologists,  one      las  received 
much  a     .on  is  that  of  heat  diffusion  in  the  a     ;ere.  *'tuch  of  this 
effort  has  bee   -  rected  towar d         •   Lurnal  temperature  wave,  and  in 
every  attempt,  the  concept  of  eddy  diffusivity  has  been  retained.  The 
classical  theory  of  daily  temperature  variations  has  bee::.  liscussed  by 
Sutton  (7),  while  the  more  recent  contributions  have  been      ,ri2ed 
by  Staley  (6) . 

■  Taylor  heat -diffusion  equation,  which  states  that  the 
turbulent  flux  of  heat  is  proportional  to  the  gradient  of  potential 


i   iture,  may  be  '.written: 


K 


Here  f  represents  the  time;  z,   the  height;  and  ^7  (z,t),  the  potential 
tperature  deviation  From  a  mean  value.   Generally,  the  coefficient  of 
edd;}  ■  ■'  'fusivity,  '-.'.,    is  a  function  of  both  height  and  time  in  an;y 
pa  rt  L c  ular  locat ion . 

In  many  of  the  past  studies,  K  has  been  assumed  to  vary,  (a) 
linearly  wit!  hei  ;ht,  i '.  )  as  a  power  of  height,  and  (0)  as  a  bounde 
exponenl  Lai  function  of  height.   In  each  of  these  case;       i  tun 
functional  form  of  K  was  assumed  to  hold  throu  There. 

Recently,  however,     ries  (l)  has  studied 

case  of  en    atmosphere  consisting  of  several  layers, 

in  each  of  which  K  is  expressed  as  a  different  function 
03"  height. 

For  an  example  he  proposed  three  layers;  a  laminar  sublayer,  e 

turbulent  boundrj    r    in      ree"  atmospheric  layer.   In      il, 

each  of  these  layers,       par   ilar,  the  turbulent  boundary 

Layer,  can  be  further  subdivided  Lnto  layers,   he  also  gave 
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(1) 


analytical  solutions  for  several  particular  cases;  (a)  K  ■  constant, 
(b)  K  =  <Xfz+z  )  ,  and  (c)  K  =  a[l  -  b  exp(-cz2J  .  Case  (a)  yields  a 
solution  In  terras  of  a  rather  complicated  exponential  function,  (b), 
a  solution  in  terms  of  Vessel  functions  of  the  first  kind,  and  (c), 
a  solution  in  terms  of  a  hyper^oraetric  function.  However,  no 
numerical  values  were  given,  and  the  results  are  not  suitable  for 
practical  purposes. 

It  is  the  purpose  of  this  investigation  to  present  a  numerical 
solution  of  2qn.  (l),  with  appropriate  boundary  conditions,  for  a 
Darticular  case  of  an  atmosphere  composed  of  an  arbitrary  number  of 
layers . 


2.   The  Heat  Diffusivity  Coefficient. 

From  physical  considerations,  K,  the  heat  diffusivity,  may  be 
expected  to  increase  in  the  lover  layers,  but  should  not  continue  to 
increase  indefinitely  with  height.  Hence,  a  bounded  diffusivity 
coefficient,  which  increases  with  height  from  the  surface  to  some  fixed 
level  and  then  decreases  with  height  to  some  residual  value  at  higher 
levels,  would  seem  more  suitable  than  either  the  case  of  K  increasing 
linearly  with  height  or  of  K  Increasing  as  a  simple  power  of  height. 
Upon  examination  of  Fig.  (1),  which  Is  based  on  observations  taken 
at  the  micrometeorological  tower  at  the  University  of  Washington  (2), 
(3)  and  Mildner's  pilot  balloon  observations  taken  near  Leipzig  (5), 
it  is  seen  that  no  single  functional  form  of  K  will  suffice.  The 
values  of  K  from  the  University  of  Washington  data  were  obtained 
by  methods  based  on  energy  continuity  at  the  earth's  surface, 
while  those  from  the  Leipzig  data  were  actually  values  of  eddy 
viscosity.  These  values  of  eddy  viscosity  were  used  as  an 
approximation  to  eddy  diffusivity,  since  they  showed  the  desired 
variation  with  height,  and  since  values  of  eddy  diffusivity  at 
sufficiently  great  heights  are  difficult  to  obtain.  Upon 
dividing  the  atmosphere  Into  four  layers,  and  then  fitting  functions 
of  height  to  each  layer,  the  following  form  of  K  resulted: {see  Fig. 

^G,  Me***.    (2) 

The  quantities  a^,a^,a_,a  ,a~,a,,  and  h  are  constants,  appropriate 
values  of  which  will  be  assigned  later.  The  discontinuity  which 
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occurs  at  z  =  228  meters  is  resolved  by  defining ^^T  at  this  level  as  zero, 

Haltiner  (k)   suggests  that  K   should  have  a  'liurnal  variation,  and  that 

it  appears  reasonable  to  expect  that  K  would  more  or  less  follow 
the  surface  tempo rature  wave,  reaching  maximum  and  minimum  values 
near  the  time  when  the  temperature  is  maximum  and  minimum, 
respectively. 

The  eddy  diffusivity,  K,  can  now  be  expressed  as  a  function  of  height  times 

a  function  of  time,  i.  e.,  K  ■  f(z)g(t),  where  g(t)  is  defined: 

g(t)  =  1  +  b(sinan  +  c  sin  2Wt)  (3) 

Here  b  and  c  are  constants,  and  CO  «=  2Tf/Q6thQQ   sec 

Mpropriate  boundry  conditions  associated  with  Eqn.  (l)  will  now  be 

chosen.  These  are  frequently  taken  to  be: 

&  =  0  f or  z  =oOy   and  all  tj  (h) 

Q   m  90  (t )  for  z  »  0;  (5 ) 

where  0Q[X)   is  defined  by: 

Qit)   -  d'sinO/t  +  c  sin  2*i>t)  (6) 

c 

Here  d  16  a  constant  and  c  and CU  are  the  same  as  in  Eqn.  (j). 


3.   Transformation  of  Coordinates. 

Tn  order  to  olace  the  basic  equations  into  a  form  which  will  give 
a  broader  interpretation  of  the  results,  and  to  scale  these  equations 
for  corrmutational  purposes,  let: 

z  =  lOOq  O  ,   and  t  =  JoOO't  .  (7 ) 

Here  q  is  a  constant,  while  ^  and  c^  are  the  new  variables.  When  t  is 
in  seconds,  the  units  of '<-  are  hours,  and  when  z  is  in  centimeters,  the 
units  of  cr  are  in  "q-meters";  i.  e.,  in  1-neter  units  when  q  =  1,  2- 
raeter  units  when  q  ■  2,  etc.  With  this  transformation,  (7),  -qs.  (l) 
and  (2)  become: 


h   =ffa-g(t)*   fit) 


-?  (9) 


<%i<r-    -f-  a 


Here  4^  must  now  be  taken  as  2~7T/2U,  and; 

a'  =  a,  (lOOq)2,         a  V  =  az flOOq), 
aj.  =  a^  (lOOq),  h*  »  h/lOOq. 

The  boundary  conditions,  Eqs.  (h)   and  (5),  remain  identical  in  form. 

Both  observation  and  theory  indicate  that  the  most  pronounced 
variations  of  temperature  occur  near  the  surface  of  the  earth. 
Therefore,  a  small  grid  distance  is  desirable  at  low  levels,  while 
at  higher  levels,  a  small  grid  distance  is  not  necessarily  needed. 
This  suggests  that  the  following  transformation  of  the  vertical 
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(8) 


coordinate  may  be  useful: 

s  a  lnCT.  (10) 

With  this  transformation,  the  lower  boundary  condition,  Eqn.  (5),  may 
apply  at  a  height  of  q  meters  corresponding  to  <=r  =  1,  however,  to  reduce 
surface  effects,  it  is  convenient  to  select  <3~  =  2.  Actually  the  vertical 
coordinate  system  is  arbitrary,  and  can  be  chosen  to  have  its  base  at  any 
level.  By  utilizing  Eqn.  (10),  Eqs.  (8)  and  (9)  now  become: 

j 


iit 


e*[K& +(&-*)!£] 


K^     to^fr)  =&L(t) 


U/ezS 

q'  e'-    *■< 


i 


Here  a y   is  taken  as  10e,  and  the  other  constants,  a'  ,  a'  ,  a  ,  a  , 
and  a*,  must  be  scaled  by  a  factor  of  10"  . 


(12) 
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k.       Finite  Difference  Equations. 

In  order  to  obtain  a  numerical  solution,  the  derivatives  of  (P 
with  respect  to  s  and  £  in  Eqn.  (ll)  must  be  replaced  by  appropriate 
finite  difference  forms.  The  problem  is  then  reduced  to  that  of 
solving  a  system  of  linear  algebraic  equations  in  the  values  of  (9 
over  a  grid  of  points  covering  the  desired  range.  Expressing  ~j  ^~a- 
and  jg  in  Eqn.  (ll)  as  finite  difference  ratios,  we  obtain: 

d  <9C    (9.  5i  i7    .  is  f    /  £  _z*>,  +  ^ ..  \ 


i-m-toX^itT*")] 


Here  the  subscript  1  designates  the  i 'th  level  of  the  vertical  grid 
at  height  i  j?  ,  where  $  is  the  vertical  distance  between  grid  points  in 
units  of  s,  that  is,  a  pure  number,  since  s  is  dimensionless . 

The  partial  differential  equation,  Eqn.  (l),  has  now  been  reduced 
to  a  system  of  ordinary  differential  equations  which  will  be 
integrated  numerically  by  the  Runge-Kutta-Gill  method. 


(13) 


5.   The  Constants. 

Text,  appropriate  valuer,  of  J( ,   n,  and  A  t,   as  well  as  the  constants 
a/  *   aJ»  etc,»  wi  11-  now  be  chosen.   Since  the  vertical  coordinate  Is 
s  =  In  a ",  .9  -  A   s,  a  constant  .  With  a  choice  of  X -   In  2,  the  i'tn 
equation  of  the  system  defined  bv  Eqn .  (15 )  applies  at  the  height 
2  q-meters.  Therefore,  the  successive  values  of  i,  beginning  with 
the  lower  boundary  condition,  anply  to  the  levels,  2,  U,  8,  16,  32, 6U, 
etc.  q-meters.  To  reduce  the  integration  time,  yet  still  give  a  fairly 
detailed  nicture  of  the  vertical  structure,  the  upper  boundary  condition, 
r-'qn .  (h)f   is  chosen  to  apply  at  the  height  of  102^  q-meters .   For  this 
choice,  <^r  ~^;   and  the  system,  Sqn  (l~>),   to  be  solved  consists  of 
eight  simultaneous  ordinary  differential  equations,  together  with 
'■  qi: .    (6)  for  &0  .     A  numerical  solution  for  tins  system  was  then 
obtained  by  the  Runge-Kutta-Gill  method  on  a  National  Cash  Register 
102A  electronic  computer.  The  choice  of  an  appropriate  time  interval 
lepended  intimately  upon  the  size  of  the  eddy  diffusivity.  As  the 
latter  increased, A  c  had  to  be  decreased  in  order  to  maintain 
computational  stability  and  reasonable  accuracy. 

The  following  numerical  values  of  the  other  constants  were 
chosen  as  being  representative: 

a'  =  0.02  (sec  l),         a^=  2050  (cm/sec), 

a j  =»  -52,500  (cm2 /sec),     aY  =  '+15,000  (cm2 /sec), 

aj  =  0.00B5,  a6=  53,000  (cm2 /sec), 

b    =    0.3,  C    =  0.5, 

d  =  0.05  (°C),  u»  =  223.  (lh) 

The  choice  of  a}  ,  n  1  ,  a-.,  a,^,  and  a/  gives  a  328-fold  increase 
of  K  with  height  from  the  surface  to  223  q-meters,  and  a  l4-fold 
iecrease  of  K  with  heigrt  above  228  q-mcter*- .  Also,  the  cnoice  of 
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b  and  c  gives  rise  to  a  21-fold  diurnal  variation  of  K  with  time.  The 
value  of  d  was  chosen  because  it  keeps  the  magnitude  of$  <0.1. 

With  the  value  of  a^  =10  ,  the  integration  interval  vas 
prohibitively  short  for  the  computer  used;  therefore  a  smaller 
value  of  a. j  =  6.25  *  104  was  used.  This  permitted  a  time  interval 
of  l'30  hour. 

The  value  of  a «  =  6.25  x  104  gives  the  following  range  of  the 
coefficient  of  eddy  diffusivity,  Eqn .  (12),  in  units  of  cm2 /sec: 

at  surface,     K /nnx~   0,15  x  lo2(C>  K/ji/a/    ~   °«'71^  q2, 

at  20  q-m,      Kmk  -  10 V,     *#,„        7-H*  q2, 

For  q  £  1,  the  values  of  K  are  somewhat  small  for  typical  atmospheric 
conditions,  especially  in  the  lowest  layers.  This  is  reflected  in 
the  results  by  a  large  amplitude  reduction  and  phase  lag  with 
increasing  height.  However,  such  values  of  K  are  not  unrealistic 
for  winter. 


at  228  q-n,     Km/   =  0.U92  x  104q2,    KmfA/   =  230  q2, 
at  14-60  q-n,     K„/jy=  0.35  x  103q2,     K  v  =  17  q2. 
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Results. 

T'hc   results  of  the  integration  are  given  by  Fir;,  (j)  and  Table  (l), 
Fig.  (5)  being  the  .graphical  interpretation  of  Table  (l).   Ii  Pi  5.  (3), 
anrolitude  reduction  and  phnse  lag,  in  hours,  are  plotted  against  height, 
in  q-meters.  As  car.  be  seen  from    .    (3)?  the  amplitude  reduction  and 
phase  lag  are  quite  large,  especially  \rher.   compared  with  observed  values. 
(See  Fig.  (*0  which  is  baser  on  observations  taken  at  the  mierometeorological 
tower  at  the  University  of  Washington).  The  large  amplitude  reduction 
an/  phase  lag  from  the  theoretical  results,  Fig.  (3),  are  due  largaly 
to  the  reduction  of  a  from  106  to  6.25  *  104.  However,  by  nroper 
selection  of  the  value  of  q,  fair  agreement  between  theoretical  and 
observed  values  can  be  obtained.  Fig.  (3)  also  shows  the  amplitude 
reduction  and  phase  lag  for  the  minimum  at  three  levels.  From  these 
three  levels,  it  con  be  seen  that  the  amplitude  reduction  is  less 
for  the  maximum  than  for  the  minimum,  the  latter  corresponding  to  the 
lower  nocturnal  diffusivity.  This  difference  is  appreciable  in  that  at 
the  eight  q-meter  level,  the  Tjaximum  is  about  k5°/o   of  the  surface  value 
O"""  £■>  while  the  minimum  is  only  19u/o  of  the  surface  minimum.  At  lower 
levels,  however,  this  difference  is  much  smaller.  The  phase  lag  increases 
with  heigh"1,  for  both  the  maximum  and  the  minimum,  however,  the  lag 
of  the  maximum  is  greater  than  that  of  the  minimum  at  any  particular 
level.  r:Tr.r  does  no+  seem  consistent  with  greater  amplitude  reduction  of 
minima  eommred  to  maxima.   Possibly,  these  differences  would  vary 
in  succeed  in?  maxima  and  minima,  as  greater  time  elapses  from  the 
initial  conditions. 

By  a  proper  choice  of  q,  some  comparisons  between  observed  and 
theoretical  results  may  be  made.  Table  (2)  3hows  amplitude  reduction 
and  nhase  lag  from  observations  ta:;eu  at  Ieafield,  along  with 
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theoretical  values  for  q  ■  10  and  q  ^  20.  For  both  values  of  q,  there 
is  fair  agreement  with  observed  values,  however,  both  the  theoretical 
amplitude  axid   phase  lag  increase  with  height  much  more  rapidly  than  do 
the  observed  values.  Table  (3)  compares  theory  with  observed  data  at 
the  Eiffel  Tower.  Here  q  is  taicen  to  be  60  and  33,  the  value  of  60 
comparing  to  the  summer  data  and  the  33,  to  the  winter  data.  Table  (U) 
shovs  some  observations  taken  at  the  University  of  Washington.  By  the 
choice  of  q  =  8,  the  amplitude  compares  fairly  well;  but  where  the  observed 
phase  lags  are  the  order  of  minutes,  those  from  theory  are  the  order 
of  hours.  In  Table  (5),  the  data  from  Porton  i3  compared  to  theory.  For 
the  choice  of  q  =  2.5,  the  comparison  is  fair,  but  once  again,  shows 
the  theoretical  phase  lag  and  amplitude  reduction  much  greater  than 
the  observed.  Table  (6)  compares  the  amplitude  reduction  in  summer 
and  winter  at  Ieafield  with  theory.  Here  q  =  5  was  chosen  for  winter, 
and  q  =  12,  for  summer.  For  the  winter  case  the  theoretical  and 
observed  values  compare  very  well,  but  for  summer,  the  theoretical 
amplitude  reduction  is  significantly  greater  than  the  observed  value. 
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Table. k. 
Amplitude  reduction  an       lag  for  the  case  of  q  =  8  versus 
observations  taken  at  the  University  of  Washington. 
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Table. 5. 
Amplitude  reduction  and  phase  lag  for  the  case  of  q  =  2.5  versus 
observation!  taken  at  ^orton,  England. 
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Amplitude  reduction  froj    e  leafield  data  versus  theoretical 
values  for  q  =  3  and  i  =  12. 
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7 .   Conclusions . 

A   numerical  solution  has  been  presented  for  the  diurnal  temperature 
variation  with  a  coefficient  of  eddy  diffusivity  which  is  a  function 
of  height  and  time.  By  suitable  selection  of  several  parameters, 
reasonably  good  agreement  has  been  obtained  between  theory  and 
observation. 

Improvement  in  the  results  may  possibly  be  achieved  by  improving 
the  functional  form  of  the  diffusivity  coefficient,  especially  in  the 
lowest  layers.  It  would  also  be  desirable  to  include  such  parameters  as 
surface  roughness,  stability,  wind  velocity,  etc.  With  a  diffusion 
coefficient  in  terms  of  these  parameters,  as  well  as  height  and  time, 
a  solution  for  the  diurnal  temperature  wave  would  have  a  much  broader 
application. 
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